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Summary 

The application of spectral methods, using a Chebyshev collocation scheme, to solve hydrodynamic 
stability problems is demonstrated on the Benard problem. Implementation of the Chebyshev collocation 
formulation is described. The performance of the spectral scheme is compared with that of a 2 nd order 
finite difference scheme. An exact solution to the Marangoni-Benard problem is used to evaluate the 
performance of both schemes. The error of the spectral scheme is at least seven orders of magnitude 
smaller than finite difference error for a grid resolution of N = 15 (number of points used). The 
performance of the spectral formulation far exceeded the performance of the finite difference formulation 
for this problem. The spectral scheme required only slightly more effort to set up than the 2 nd order finite 
difference scheme. This suggests that the spectral scheme may actually be faster to implement than higher 
order finite difference schemes. 




1.0 Introduction 

The theory of hydrodynamic stability has helped to explain and predict a variety of fluid 
flow phenomena. Recently it is being used to guide the modem computational fluid 
dynamicist in choosing the appropriate parameter values which are needed to simulate 
fluid flow behavior of interest (NASA TM-4569, 1994). Many current applications of 
hydrodynamic stability theory are possible because the field has benefitted greatly from 
the development and refinement of computational tools in addition to the existence of 
increasingly powerful computers. Spectral methods is one such set of tools that has been 
successfully applied to obtain high accuracy hydrodynamic stability results to previously 
intractable problems. 

The purpose of this paper is to show, by example, the use of a spectral collocation 
formulation to solve hydrodynamic stability problems. Our discussion will be confined to 
the linear stability analysis which is the foundation of hydrodynamic stability theory (Lin, 
1945). The linear stability problem ultimately reduces to a matrix eigenvalue problem, 
and the peril of the eigenvalue problem is that it requires 0(N^) operations to obtain the 
eigenvalues where the matrix is N x N. As shown herein, the high accuracy of spectral 
methods results in small N, therefore considerably less CPU time is required to solve for 
the eigenvalues when compared to finite difference methods. 
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The Benard problem is used to illustrate the implementation and performance of the 
spectral scheme. The problem is also solved using a 2nd order finite difference scheme 
which required slightly less time to implement. Results of the two numerical schemes are 
compared to the exact solution of the Marangoni-Benard problem (Pearson, 1959). 

The Benard problem is described, and the governing equations and boundary conditions 
are developed in the following section. After a brief description of the finite difference 
* scheme, the spectral collocation formulation is discussed. Results from both numerical 
schemes are then compared to an exact solution of the Marangoni-Benard problem. The 
spectral scheme yields results with considerably better accuracy using an order of 
magnitude less points than the finite difference scheme. 

2.0 Description of Benard Problem & Development of Equations 

A temperature difference is imposed normal to the free surface of a thin liquid layer of 
fluid of infinite horizontal extent and finite thickness, d, as shown in Figure 1 .The initial 
steady state or base state of the system is one of no fluid motion, with a linear 
temperature profile across the layer. The velocity and temperature profiles illustrated in 
Figure 1 can immediately be expressed as, U b = 0 and T b * = T b * 0 -pz*. Using the notation 
of Pearson (1958) and Chandrasehkar (1981), U b and T b * are respectively, the base flow 
velocity and temperature. The temperature gradient of the base state, P is defined as 
p = -dT b */dz'or p = AT d */d where AT d * = T b * 0 -T b * s . The asterisk, denotes dimensional 
quantities. The lower surface is rigid and is held at a constant temperature. The upper 
surface is free and exchanges heat with the environment . The free surface is assumed flat 
which is physically justified for many terrestrial problems. We first give the 
nondimensional form of the governing equations and in the next section we linearize 
about the base state just described in order to determine whether small disturbances to the 
base state will grow or decay. Specifically we are interested in the critical values of the 
nondimensional parameters where the change of stability occurs. 


z** 


Q* 



Figure 1 Base State For Thin Liquid Layer Of Infinite Extent 


Nondimensional forms of mass, momentum, and energy equations for an incompressible 
fluid with the Boussinesq approximation are given in equations (1) through (3). The 
derivation of these equations with the Boussinesq approximation and constant viscosity 
and their subsequent nondimensionalization are well known and we refer the interested 
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reader to Chandrasekar (1981), and Drazin & Reid ( 1 982) for details. All 
thermophysical properties are assumed constant apart from density and surface tension. 


V-U = 0 (1) 

^ = -VP + kRaPr(T-T b0 ) + PrV 2 U (2) 

Dt 

— = V 2 T (3) 

Dt 

0, T, P, t are the velocity vector, temperature, pressure, and time respectively. The 
reference values used to nondimensionalize the variable; length, velocity, temperature, 
pressure, and time are d, k 0 /d, pd, p^/d 2 , d 2 /^ * respectively, pj is the fluid 
density and k 0 is the fluid thermal diffusivity. The subscript 0 indicates that the 
properties are chosen at the lower surface temperature, T^o. The characteristic value of 
the dynamic viscosity of the fluid, p, is denoted as pj,. These reference values are 
consistent with those used in the buoyancy instability studies presented in Chandrasekhar 
(1981) and Drazin and Reid (1982), and the surface tension instability investigations of 
Pearson (1958) and Scriven and Stemling (1964). Two dimensionless groups appear in 
the momentum equation, the Prandtl number, Pr, and the Rayleigh number, Ra, which 
are defined as follows: 

Ra = *&. 

Po*0 K oPo 

is the volumetric thermal expansion coefficient and g z is gravitational acceleration in 
the negative z-direction. The dot product of the unit vector in the z direction, k, and the 
buoyancy (RaPr) term in equation (2) indicates that buoyancy only acts in the vertical 
direction. Therefore the Rayleigh number only occurs in the z-momentum equation. 

The nondimensional boundary conditions are given by equations (4) and (5). Equations 
(4a,b,c), represent the no-slip conditions and impenetrable wall condition at z=0. 

Equation (4d) is the constant temperature condition along the wall. The normal stress 
boundary condition reduces to (5a) when the free surface at z=l is assumed to be flat. 
Boundary condition (5b) is the heat flux balance at the free surface, where Q* is the 
dimensional surface heat flux to the environment and k b is the fluid thermal conductivity. 
Equations (5c and 5d) are the tangential force balances along the free surface, in the x and 
y directions, respectively. 


At z = 0; 
At z = 1 ; 


U(0) = (U x ,U,,U z ) = 0, T(0) = T, 

u - (i)=0; S + S'°' 


bO 


( au, | au, \ | ( au, | au, 

<3x dz J y dy dz 


P 


= MaV„ T 


(4a,b,c,d) 
(5a, b) 

(5c, d) 
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The operator, V„ , is the surface gradient defined as i — + j — where, i and j are unit 

dx dy 


vectors in the x and y directions respectively. The Marangoni number, Ma, which occurs 

in equations (5c,d) is defined as: Ma = ^ . The parameter, y *, is defined as 

K 0 p 0 ST* 

and is often referred to as the temperature variation of surface tension (Nield 1964 and 

Adamson 1967) or differential coefficient of surface tension change with temperature 


(Scriven and Stemling 1964). The surface tension, a*, does not appear in our equations or 
boundary conditions since we have assumed a flat interface. Further discussion of the 
nondimensionalization of the free surface boundary conditions is found in Scriven and 
Stemling (1964), and Koschmeider (1993). 


The surface heat flux, Q*, has to be expressed in a form that is suitable for linearizing the 
heat flux boundary condition, equation (5). This is accomplished by expanding Q* about 
the base state surface temperature, T b ’ s . The first order expansion is given by equation 
(6). As previously noted, the base state varies only in the z-direction. Therefore, Q*(T bs ) 
can be re-expressed as equation (7), using Fourier's law. 


Q’=Q'CO + |j!r| 

Q'(T;)=k;^r| 

dz Ud 


(r -X.) 

£ 

=k;p 


( 6 ) 

( 7 ) 


Substituting equation (6) into equation (5b), using koP in place of Q’(tJ s ) and defining 


h* = 



> the heat flux boundary condition becomes: 


fI + l + Bi > (T,-Tj = 0 
dz 


( 8 ) 


h*d 


The dimensionless group, Bi s , is defined as Bi s = — — and is referred to as either the 

K 

surface Biot number (Pearson, 1958 and Nield, 1964) or the surface Nusselt number 
(Scriven and Stemling, 1964). 


We note that the three-dimensional mass, momentum, and energy equations are given in 
equations (1-3), yet the boundary conditions are only specified in the z-direction. After 
linearizing the problem and applying some vector operations, it is shown in the next 
section, that the governing equations and boundary equations in the x and y directions do 
not affect the stability of the base state. Equations (1, 2, 3, 4 , and 5a,c,d and 8) make up 
the system which we will linearize in the next section. 
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2.1 Linearization of The Governing Equations 


The dependent variables are written in terms of the following base flow and perturbation 
variables: U = u, T = T b + 0, AP = AP b + Ap 

After substituting for T b and VT b , the disturbance equations become: 


— = - Vp + k • Ra Pr 0 + Pr V 2 u 

(9) 

a 


— -u. = V 2 0 

(10) 

dt ' 



k is the unit vector in the z-direction shown in Figure 1 . The curl operator is applied 
twice to the momentum equation, equation (9), which yields equation (11). 


c(v 2 u) „ ' .... 

— - = k • Ra Pr V 2 0 + Pr V 4 u (11) 

at 

The first curl operation yields the vorticity equation and eliminates the pressure terms. 
The second curl operation decouples the momentum equations from each other. The x 
and y momentum equations become uncoupled from the z-momentum and the energy 
equations. The z-momentum and energy equations remain coupled through the buoyancy 
term in equation (11), the convective term in equation (10), and the tangential boundary 
condition (discussed below). Furthermore, the relevant stability parameters, Ma and Ra, 
do not appear in either the x or y momentum equation or their associated boundary 
conditions. Given these considerations, equation (11) reduces to the scalar equation in u z , 
equation (12). 

— = RaPrV 2 ,0 + PrV 4 u (12) 

dt 

The boundary conditions for the perturbed variables associated with equations (10 and 
12) are given by equations (13) through (14). 


At z = 0, 
At z = 1 , 


u = 0; % = 0; 6 = 0 

dz 

u z = 0; ^ + Bi s 0(l) = o 

dz 


{ 


d 2 u a 


V dz J 


= MaV 2 0 


(13a, b,c) 
(14a,b) 

(14c) 
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2.2 Normal Mode Analysis 


Since equations (10) and (12) are linear, we assume solutions for u z and <f) are of the 
form: 

, v i(a,x+a v y)+Xt , A , , s i(a,x+a v y)+Xt 

u 2 = w(z)e ' * ' ’ and 0 = <(>(z)e 

a x and ay are the dimensionless wavenumbers in the x and y directions, and A is the 
dimensionless frequency. Substituting these into equations (10) and (12 ) results in the 
following ordinary differential equations. 

>4>(z) - D 2 <)>(z) + a 2 <|>(z) - w(z) = 0 (15) 

a(d 2 w - a 2 w(z)) = Ra Pr(D 2 <J> -a 2 <|>) + Pr(D 4 w - 2a 2 D 2 w + a 4 w(z)) (16) 

Where D = — and a 2 = a 2 +a 2 . 

dz y 

The boundary conditions at z = 0 become: 

w(0) = 0, Dw(0) = 0, <(>(0) = 0. (17a,b,c) 

At z=l, the flat interface condition, heat flux condition, and tangential stress boundary 
condition are: 

w(l) = 0, D<)>(l) + Bi s <t>(l) = 0, D 2 w = -a 2 Ma<|)(l) (18a,b,c) 

Equations (15 through 18) are solved to determine whether the velocity and temperature 
disturbances grow or decay for given combinations of the relevant parameters. The 
relevant parameters are Ma, Ra, and a. Our problem is also referred to as a temporally 
developing flow problem since the disturbance growth or decay is in time. For temporally 
developing flows, a x and ay are real and the eigenvalue. A, is complex. If the real part of 
A. is positive the disturbance grows, if the real part of A. is negative the disturbance decays 
in time and if A, is zero, the disturbance persists unchanged in time. 

3.0 Discrete Formulations 

Two discrete formulations will be applied to the Benard problem, a 2nd order finite 
difference scheme, and a spectral collocation scheme. Irrespective of the discrete 
formulation the goal is construct a set of linear equations in form of the general 
eigenvalue problem, Ax = ABx. Once the eigenvalue problem is setup, solution 
mechanics are identical. If B is cheaply invertible, it usually pays to reduce the problem 
to a regular eigenvalue problem of the form Cx =A.x, where C = B*^A. In this study we 
used standard QR and QZ eigenvalue subroutines from the IMSL library to solve for the 
eigenvalues. 
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3.1 2nd Order Finite Difference Scheme 


Equations (15 through 18) were discretized using a standard central difference scheme. 
Boundary conditions result in solving 2N-1 equations for this formulation. The 
discretized governing equations, equations (19) and (20), were arranged in the form 
Az = XBz, which is the generalized eigenvalue problem. Coefficients, a through g, and r 
are constants. 


aWj., +bw j _, +cwj +dw itl +ew i+2 + r<|) j = X(fw M +gw, +fw i+1 ) (19) 

f<t>i-i+g<t>i +f< J>i + i+Wi =^(<10 (20) 


Here the boundary conditions for Eq (15) are applied only to the i = 0, 1, N-l, N 
equations and for Eq (16) only to i = 0 and i = N. 

B is a nonsingular matrix so it is possible to reduce the system to a regular eigenvalue 
problem of the form Cz = B _I A = Xz. Assuming a flat interface ensures that B is a 
tridiagonal matrix which can efficiently be inverted using a tridiagonal solver. The 
problem was discretized in terms of one fourth order equation, and one second order 
equation, which yields A and B matrices of rank 2N+2. Three of the six boundary 
conditions are Dirichlet boundary conditions which reduce the A and B matrices to rank 
2N-1. 


3.2 Chebyshev Collocation Spectral Scheme 

The key to all spectral techniques lies in the possibility of expanding smooth functions in 
terms of rapidly converging sums of certain orthogonal basis functions. For example, 
consider any reasonable function f(x) defined in the domain -1 < x < 1 (see Canuto et. 
al. for a precise definition of "reasonable"). The function can be represented as a sum of 
Chebyshev polynomials, T n (x), of the form: 

f(x) = £cr„(x) (21) 

n*0 

The crucial thing is that the sum converges very rapidly if f(x) is smooth so one can 
truncate it at N terms and accurately represent the function with a minimal set of numbers 
n = o n}- Such an expansion can be viewed as a very efficient and only slightly lossy 

compression technique for functions. 

Pure spectral methods proceed by expanding the unknowns in terms of truncated sums of 
certain polynomials having excellent convergence properties (often simple combinations 
of Chebyshev polynomials that automatically account for any boundary conditions that 
must be satisfied by the function). The sums are then substituted into the differential 
equation and the coefficients are picked to minimize the residual. The fundamental 
quantities of interest in this procedure are the coefficients in the expansions of the 
dependent variables. 
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Spectral collocation methods on the other hand, concentrate directly on the physical space 
representation of the unknowns and as a consequence are more easily understood by the 
naive user. For example, in a collocation technique our hypothetical function, f(x) is not 
stored as {^ : n = o n}, but instead as = f( x , ) : i = 0,...,n}. The exact correspondence 

between the two representations is maintained by choosing the physical space grid {xj 
in an optimal fashion that is related to one of the Gaussian integration formulas. A typical 
formula of choice for Chebyshev expansions on the domain [-1,1] is the Gauss-Labatto 
grid, x, = cos(m/N). The spectral space and physical space representations can be 
interchanged with essentially no error (except perhaps for aliasing errors). Moreover, if 
the expansion is in terms of Chebyshev polynomials or trigonometric functions, the 
transformations to and from spectral space can be carried out efficiently by using Fourier 
Transforms (FFT's). 

Solving differential equations obviously requires that the derivatives are evaluated. One 
method of evaluating {(-s f(x,)} is to proceed as 

That is, one first transforms to spectral space where a derivative is taken rapidly by using 
some simple properties of the basis functions. The new series produced in this fashion is 
then transformed back to a physical space representation. In the case of trigonometric or 
Chebyshev expansions, the procedure is dominated by the FFT's used in the 
transformations, so the total cost is O(NlogN) operations. 

There is a mathematically equivalent approach which uses matrix-vector multiplies to 
express 

j-0 

where the elements of the derivative matrix D can be found in spectral texts such as 
Canuto et. al. To evaluate the derivatives on the entire grid using this method will take 
0(N 2 ) operations. However, the matrix-vector multiply approach is the only one possible 

for eigenvalue problems where the aim is to turn the linear differential operator into the 
equivalent matrix operator on the discrete grid. Thus for example, the continuous 
equation, Xf = f\ becomes the discrete equation, Xf = D 2 f, so that in theory one simply 
fills and then squares the D matrix before feeding it to a standard matrix eigenvalue . In 
practice, the greatest programming labor is in the implementation the boundary 
conditions. 

For the Benard problem, it is convenient to define the matrix operator 

L = D 2 -a 2 I (23) 
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where I is the identity matrix. The continuous equations, (15) and (16) then become the 
discrete equations, (24) and (25) 

Pr(L 2 w-a 2 RaI<|) J ) = ?Xw (24) 

L<|> + Iw = ?J<t> (25) 

and these equations can be recast into the standard form Au = XBu given by equation 
(26). 


PrL 2 

-orRaPrl 

w 

= x 

~L 

0 w 

I 

L 

A 


_0 

iJlA 


(26) 


This translation of a continuous problem into a discrete one is very natural and can be 
carried out even more rapidly than the corresponding process for a finite difference 
scheme. However, this matrix eigenvalue problem as it stands does not take the boundary 
conditions into account. Most of the coding complexity that is present in spectral 
techniques (which by nature are global approximations) arises because of the need to 
implement boundary data (which are local point conditions). 

The boundary conditions and governing equation are first mapped from the z variable in 
the domain [0,1] to x defined on [-1,1] by x = 2z-l. The mapped boundary conditions 
become: 


At x = -1 (i=n): w n = £d n ,w, =<(>,= 0 (27) 

i=0 

N N 

At x = +1 (i=0): w 0 = Bi s <|> 0 +XDoi<l>i = ct 2 Ma<|> 0 £D 2 w, =0 (28) 

i =0 1=0 

We immediately see that w 0 = w N = <(> N = 0. The remaining equations can be used to 
simultaneously solve for w,, w N _,, and <j> 0 in terms of the other w/s and 4>,’s as shown 
in equation (29). Equations (29) reveal that elements, w, and w N _,, are coupled to the 
<(>. ’ s through the last boundary condition in equation (28) while <j> 0 remains uncoupled 
from the w/s. 


N-2 N-2 N-l 

W l =Z C ' W ' +< f > 0’ W N-1 =X d i W 1 +< t’0> <(>0=2]^! ( 29 ) 

1=2 i=l 

The boundary condition information is used to reduce the rank of the eigenvalue problem 
which is given by equation (30). 


2N+1 


2N+I 


j=0 j=0 


(30) 
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The first N+l components of the vector u here corresponds to {w ( :i = 0,...,N} and the 
last N+l components corresponds to {^.’i = 0,...,N}. Applying the six boundary 
conditions, we can eliminate the rows corresponding to w 0 , w, , w N _, , w N , <(> 0 , and <j) N and 
expand equation (30) as given below. 


2N+1 


2>u u i = 


J-2 


j*N-l.N.N+l 



2N+1 

BjjUi + B iN _ 1 u N _, + B iN+1 u N+) + 

Z B .i u , 


j-2 


(31) 


j*N-1,N,N+I J 


where we have already used the data. u 0 = u N = u, N+ , = 0 (which corresponds to 

w 0 = w N = <|) N = 0). Using the remaining boundary conditions, equations (29), take the 

form 


N-2 2N N^2 2N_ 

u . = Z c J u j= Zvv u N-. = Z d j u j= 2- d j u j’ 


J*2 


N-l 


J=2 

j*N-l,N,N+l 
2N 


j*2 


J'2 

j*N-I.N,N+l 


u N+l =Ze jUj = X e j u j 


j*2 

j*N-l,N.N + l 


d: 


where C) = \ 1 & d J= ] J & ^ = 


fo 


'j-N-3 


'j-N-3 


'j-N-3 


The matrix eigenvalue problem can be rewritten as 


ifj = 2,...,N-2 
ifj = N + 2,...,2N 


(32) 


(33) 


2N+I 2N+1 

Za, U j =X Z BijUj (34) 

j=l j*l 

j#N-l,N,N+l j^N-l,N,N+l 

where A {J = A y +A il c j +A iN .,d J +A jN+1 e j and B u = B (j + B iI c J + 8^.,^ + 3^,^ . The 

global matrix problem is finally transformed to a reduced eigenvalue problem Au = aBu 
where the matrices are (2N -4) x (2N -4) and all six boundary conditions have been 
incorporated into the problem. 

The principle difficulty in using spectral collocation techniques for solving stability 
problems is the implementation of derivative boundary conditions. A condition such as 
w( 1 ) = 0 is not a problem as all that is required is the reduction of the global matrices by 
eliminating one row and one column. Derivative data on the other hand, results in altering 
all the elements of the matrices due to the global nature of the underlying series 
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approximation. All data contributes to the value of the derivative at each point 
Nevertheless, for the Benard problem, the spectral collocation technique is only slightly 
more difficult to code than a second order finite difference scheme. As is observed in the 
subsequent results, the additional coding effort is amply rewarded. 

4.0 Results 

Under certain conditions, exact solutions to equations (15) through (18) have been 
obtained. Pearson (1958) derived an exact solution to the Marangoni-Benard problem 
(Ra = 0) for the case of neutral stability, X = 0. His solution reduces to equation (35) for 
an insulated free surface, Bi = 0. The critical value of the Marangoni number, equation 
(35), versus the wavenumber is shown in Fig. 3 and is referred to as a neutral stability 
curve since X = 0 for all points along the curve. For values of Ma above this curve are 
unstable since infinitesimal disturbances. Our objective is to use the above exact result to 
investigate the accuracy of the aforementioned discrete formulations, so we do not 
consider alternative exact or approximate solutions which exist for the general problem. 
The physical interpretation of these results in addition to results from other exact or 
approximate solutions to the Benard problem are discussed in Pearson 1958, Scriven & 
Stemling 1964, Smith 1966, Chandrasekhar 1981, and most of the other references cited 
in section 6.0. We now compare the numerical results to the exact solution. 


Ma c = 


8a 2 cosh(a)(a - sinh(a) cosh(a)) 
a 3 cosh(a) - sinh 3 (a) 


(35) 



Figure 2 Marangoni-Benard Neutral Stability Curve, 
Exact Solution for Bi=0 (Pearson, 1958) 
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Neutral stability curves for the Marangoni-Benard problem which were generated using a 
2nd order finite difference scheme and spectral scheme are shown in Figs. 3 and 4, 
respectively. The number of points across the fluid depth, N, (in the z-direction) 
represents the spatial resolution used to generate a given curve. In Fig. 3 the neutral 
stability curve converges to the exact solution as the spatial resolution increases from 
N = 4 to N = 100. The N = 50 and N = 100 curves are visually indistinguishable from the 
exact solution. Fig. 4 reveals that the neutral stability curves computed using the spectral 
* formulation also converge to the exact solution as the spatial resolution increases. The 
spectrally generated neutral stability curves shown in Fig. 4 are visually identical to the 
exact solution for spatial resolutions as low as N = 10. In both Figs. 3 and 4, the 
numerically generated neutral stability curves tend to diverge from the exact solution with 
increasing wavenumber. It is also observed that the finite difference solution converges 
from above the exact solution while the spectral solution converges from below exact 
neutral stability curve. 



Figure 3 Marangoni-Benard Neutral Stability Curves 
Computed Using A 2nd Order Finite Difference Scheme. 
N Is The Number Points Through The Fluid Layer. 



Figure 4 Marangoni Neutral Stability Curves 
Computed Using A Spectral Collocation Scheme. N 
The Number Points Through The Fluid Layer. 


The error in the Marangoni number for the finite difference and spectral schemes is 
plotted as a function of wavenumber in Figs. 5 and 6, respectively. In both figures, error 
is plotted using a logarithmic scale while the wavenumber, a, is plotted with a linear 
scale on the abscissa. The error (ordinate) range differs between the two figures so that 
the error characteristics of each discrete scheme could be observed. The error is defined 
as Ma N - Ma CXJC1 where Majq is the Marangoni number computed from a discrete 
Ma eua 

formulation for a given spatial resolution (N points) and Mae Xact is computed using 
equation (35). Both discrete schemes are observed in Figs. 1 through 4 to converge to the 
exact solution as the spatial resolution, N, increases. The finite difference errors for each 
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curve increase approximately an order of magnitude over the given wavenumber range as 
observed in Figure 5, while the spectral error shown in Figure 6 increases four to five 
orders of magnitude with wavenumber. For N greater than approximately seven, the 
spectral error remains considerably less than the finite difference error for the range of a 
considered. 

N = 4 

N = 5 

N = 7 

N = 10 

N = 15 

N = 50 

N = 100 

0 1 2 3 4 5 6 7 

a 

Figure 5 Marangoni Number Error vs. Wavenumber For The 2nd Order Finite Difference Solution. 
N Is The Number Of Points Through The Fluid Layer 


N = 4 
N = 5 
N = 7 
N = 10 
N = 15 

a 

Figure 6 Marangoni Number Error vs. Wavenumber For The Spectral Collocation Scheme. N Is 
The Number Of Points Through The Fluid Layer 
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Selected error values of the two schemes for wavenumbers of 2 and 5 are tabulated in 
Tables 1 and 2, respectively. The errors for , a = 2, the critical wavenumber were smaller 
than those of the larger wavenumber, a = 5. In Fig. 7, the Ma error is shown on a log-log 
plot as function of the spatial resolution (N grid points) for both discrete formulations . 
Comparing Ma error values at a wavenumber of 2 and N = 10, the error for the spectral 
scheme is seen to be five orders of magnitude smaller than the finite difference error. 
Furthermore, after increasing the spatial resolution of the finite difference scheme to 
, N = 100, still gives a spectral error for N = 10 that is 3 orders of magnitude smaller. The 
reduction in error for the finite difference scheme is essentially proportional to N 2 , as 
expected since the scheme is 2 nd order accurate. A slope of -2.02 was computed for the 
finite difference curve in Fig. 7 which is within 1% of the expected value of 2. The slope 
was computed from a least squares fit of the finite difference data in Table 1 . The error 
for the spectral formulation is expected to decrease exponentially with increasing N 
(Boyd, 1989, Canuto et. al., 1987). However the error results in Table 2 show that the 
exponential rate of convergence is exceeded for this particular problem. Fig. 7 vividly 
illustrates that the spectral scheme results in a significant reduction in error with 
considerably fewer grid points than the central difference scheme for this particular 
problem. The spectral formulation has also been shown to out perform finite difference 
methods when applied to other hydrodynamic stability problems (Canuto et. al., 1987, 
Boyd, 1989). The exceptional performance (greater than exponential convergence) of the 
spectral collocation scheme for this problem was not anticipated by the authors. 


Table 1 Selected Marangoni Number Errors 
For Wavenumber, a=2 


Spatial resolution 
N 

Finite 

Difference 

Spectral 

5 

1. 085X10* 1 


10 



15 

pfliuma 

4.529x10-* 1 

50 

1.025x1 0' 3 


100 

2.56 lxl 0* 4 



Table 2 Selected Marangoni Number Errors 


For Wavenumber, a=5 


Spatial resolution 
N 

Finite 

Difference 

Spectral 

5 

5.756x10-* 

2.808x10-* 

10 

1.309x10'* 

2.525x1 O' 4 

15 

5.710xl0* 2 

6.434x1 O' 9 

50 

5.075x1 0* 3 


too 

1.268x1 0' 3 



As stated throughout, the ability to reduce the size of N is crucial to the eigenvalue 
problem. Ax = a.Bx. Inverting B takes 0(N^) operations; the matrix multiplication of 
C = B*'A requires 0(N^); and solving the regular eigenvalue problem, Cx = kx, requires 
OCN^) operations. Neglecting all other operations than those identified above, for a grid 
resolution of N=10, it requires appropriately 3000 operations to compute the eigenvalues 
while it requires 0(3x1 0^) operations for N=100. The number of iterations required to 
converge to Mac at one wavenumber is 0(10), ie.. the matrix eigenvalue problem is 
solved approximately ten times for each wavenumber. 
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Figure 7 Error in Ma c (a = 2) As A Function Of The Number of Points 
Across The Fluid Layer, N 


5.0 Concluding Remarks 

A spectral scheme, the Chebyshev collocation formulation was used to perform a 
hydrodynamic (linear) stability analysis of the Benard problem. The problem reduces to a 
I generalized eigenvalue problem. Ax = XBx, which can be reduced to a regular eigenvalue 
problem, Cx = kx by inverting B. Implementation of the spectral scheme was described. 
There is a bit of a learning curve that must initially be overcome to comfortably setup the 
spectral formulation if one has no previous experience with spectral methods. Afterwards, 
the spectral scheme requires only slightly more time to set up than the 2nd order finite 
difference scheme and is likely to be easier to program than higher order finite difference 
I schemes. A comparison of the results from the spectral and finite difference scheme 
reveals that the spectral scheme out performs the finite difference scheme by a 
considerable margin. The error of the spectral scheme is at least three orders of magnitude 
smaller than the finite difference error for N = 10 and seven orders of magnitude smaller 

than finite difference error for N = 1 5. 

i 
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